pro javadatabase

	joLDB = connect_db()
	;assume success

	;490038.976	5263768.19

	centrex = 490038.976
	centrey = 5263768.19
	;half side lengths
	sizebig = 150
	sizesmall = 10
	height = 500 ;0 to height
;	boxbig = "BOX3D(" + strtrim(centrex - sizebig,1) + " " + strtrim(centrey -sizebig,1) + " 0, " + strtrim(centrex + sizebig,1) + " " + strtrim(centrey + sizebig,1) + " 500)"
;	boxsmall = "BOX3D(" + strtrim(centrex - 50,1) + " " + strtrim(centrey - 50,1) + " 0, " + strtrim(centrex + 50,1) + " " + strtrim(centrey + 50,1) + " 500)"
;	print, boxbig, boxsmall
;
;	newq = "select point_geom from points where point_geom && '" + boxbig + "'::box3d AND class = 2 AND z(point_geom) <= 500"
;	newqq = "select point_geom from points where point_geom && '" + boxsmall +"'::box3d AND class = 0"
;	print, newq
;	print, newqq
;	print, joLDB->executeQuery(newq)
;
;	ground = joLDB->getAllPoints()


	points = fetchc(centrex, centrey, 100, dbc=joLDB)

;	all = "select point_geom, extra from points, points_extra where tile_id = 'pt000166' AND points.id = points_extra.id LIMIT 50000"
;	some = "select point_geom, extra from points, points_extra where point_geom && 'BOX3D(562482.77 5405821.0 0, 562832.77 5406171.0 100)'::box3d AND z(point_geom) <= 600 AND points.id = points_extra.id"
;
;
;	x = joLDB->executeQuery(some)
;	points = joLDB->getAllPoints()
;
;	joLDB -> ResultSetReset
;	metas = joLDB->getColumnAll(2)
;
;	meta = dblarr(n_elements(metas),9)
;		v1 = 0.0d
;		v2 = 0.0d
;		v3 = 0.0d
;		v4 = 0.0d
;		v5 = 0.0d
;		v6 = 0.0d
;		v7 = 0.0d
;		v8 = 0.0d
;		v9 = 0.0d
;	for i=0L, n_elements(metas) - 1 do begin
;		reads, metas[i], v1, v2, v3, v4, v5, v6, v7, v8, v9
;		meta[i,0] = v1
;		meta[i,1] = v2
;		meta[i,2] = v3
;		meta[i,3] = v4
;		meta[i,4] = v5
;		meta[i,5] = v6
;		meta[i,6] = v7
;		meta[i,7] = v8
;		meta[i,8] = v9
;	endfor
;
;	save, points, meta, filename="D:\sjwelch\pntmeta.dat"

	restore, "D:\sjwelch\pntmeta.dat"

	;colour table stuff
	loadct, 13
	tvlct, r, g, b, /get
	rgbtable=BINDGEN(3,256)
 	rgbtable[0,*]=r
    rgbtable[1,*]=g
    rgbtable[2,*]=b

	vertc = intarr(n_elements(points[*,0]))
;	maxh = max(points[*,2])
;	minh = min(points[*,2])
;	vertc = ((points[*,2] - minh) / (maxh - minh)) * 255
	f = 4
	maxh = max(meta[*,f])
	minh = min(meta[*,f])
	vertc = ((meta[*,f] - minh) / (maxh - minh)) * 255

	ground = points[where(meta[*,4] eq 2),*]

	iplot, transpose(ground), /scatter, identifier=tool_id, /no_savepromt;, rgb_table=rgbtable, vert_colors=vertc
	;isurface, array[0,*], array[1,*], array[2,*], overplot=tool_id

	STOP
	print, joLDB->executeQuery(newqq)
	nong = joLDB->getAllPoints()
	iplot, transpose(nong), /scatter, overplot=tool_id, /no_saveprompt

   	ods = currentds(tool_id)

	joLDB->close


	ground = transpose(ground)
	nong = transpose(nong)

	;BUILD DTM
	gridsize = 1024

	dtm = builddtm(ground, gridsize)
	dtmbounds = [min(ground[0,*]), max(ground[0,*]), min(ground[1,*]), max(ground[1,*]), min(ground[2,*]), max(ground[2,*])]
	chm = buildchm(nong, dtm, dtmbounds)

	save, chm, FILENAME="D:\sjwelch\chm.dat"

	mesh_obj, 0, verts, polys, ground
	vism = obj_new('IDLitVisPolygon')
	vism-> setproperty, data=ground, polygons=polys
	vism -> setproperty, linestyle=6, fill_background=1, fill_color=[255,0,0], color=[0,0,0], noClose=0, shading=0
	ods->add, vism

	palette = transpose(rebin(bindgen(256), 256, 3))
	visi = obj_new('IDLitVisImage', image=dtm, visualization_palette=palette)
	xr = (max(ground[0,*]) - min(ground[0,*]))
	yr = (max(ground[1,*]) - min(ground[1,*]))
	xs = xr / gridsize
	ys = yr / gridsize
	;iimage, dtm, overplot=tool_id
	visi -> setproperty, xorigin=min(ground[0,*]), yorigin=min(ground[1,*]), zvalue=min(ground[2,*]), pixel_xsize=xs, pixel_ysize=ys
	ods->add, visi



	points = dblarr(gridsize,gridsize)
	count = dblarr(gridsize,gridsize)
	average = dblarr(gridsize,gridsize)
	sclground = [(ground[0,*] - min(ground[0,*])) / (max(ground[0,*]) - min(ground[0,*])), (ground[1,*] - min(ground[1,*])) / (max(ground[1,*]) - min(ground[1,*])), (ground[2,*] - min(ground[2,*])) / (max(ground[2,*]) - min(ground[2,*]))]
	sclnong = [(nong[0,*] - min(ground[0,*])) / (max(ground[0,*]) - min(ground[0,*])), (nong[1,*] - min(ground[1,*])) / (max(ground[1,*]) - min(ground[1,*])), (nong[2,*] - min(ground[2,*])) / (max(ground[2,*]) - min(ground[2,*]))]

	points[sclground[0,*] * gridsize - 1, sclground[1,*] * gridsize - 1] += ground[2,*]
	count[sclground[0,*] * gridsize - 1, sclground[1,*] * gridsize - 1]++
	count[where(count lt 1)] = 1
	average = points / count

	;save, dtm, sclground, ground, nong, points, average, count, FILENAME="D:\sjwelch\idlvars.sav"

	dtmorpoints = dtm OR average
	absnong = nong
	absnong[2,*] = nong[2,*] - dtmorpoints[sclnong[0,*] * gridsize - 1, sclnong[1,*] * gridsize - 1]
	;iplot, absnong, /scatter, /no_savepromt;, overplot=tool_id

	mind = [min(absnong[0,*]), min(absnong[1,*]), min(absnong[2,*])]
	maxd = [max(absnong[0,*]), max(absnong[1,*]), max(absnong[2,*])]

	;BUILD CHM

	;build volume to extract CHM (finding the highest point in each stack=neighbourhood)
	xrnong = max(nong[0,*]) - min(nong[0,*])
	yrnong = max(nong[1,*]) - min(nong[1,*])
	vol = ulonarr((xrnong/xr)*256, (yrnong/yr)*256, 512)
	;vol = build_volume(SIZE(vol, /dimensions), absnong, mind, maxd, vol=vol)
	;ivolume, vol

	volvs = size(vol, /dimensions)
	hm=ulonarr(volvs[0], volvs[1])
	arrin = array_indices(vol, where(vol gt 0))
	triangulate, arrin[0,*], arrin[1,*], arrtri
	chm = trigrid(arrin[0,*], arrin[1,*], arrin[2,*], arrtri, NX=volvs[0], NY=volvs[1])
	hm[*] = 0
	hm[arrin[0,*], arrin[1,*]] = hm[arrin[0,*], arrin[1,*]] > arrin[2,*]
	tvscl, hm

	;now build trigrid of these points
	triangulate, absnong[0,*], absnong[1,*], nongtri
	xrnong = max(nong[0,*]) - min(nong[0,*])
	yrnong = max(nong[1,*]) - min(nong[1,*])
	chm = trigrid(absnong[0,*], absnong[1,*], absnong[2,*], nongtri, NX=(xrnong/xr)*gridsize, NY=(yrnong/yr)*gridsize)
	isurface, chm
	;tvscl, smooth(chm, 3, /EDGE_TRUNCATE)

	save, chm, FILENAME="D:\sjwelch\chm.dat"


end

;calculates per vertex colours which are indices into the colour table (vals 0...255)
function calc_vert_colours, valspervert
	mnn = min(valspervert, max=mxx)

	return, ((valspervert - mnn) / (mxx - mnn)) * 255
end



function findbounds, cloud
	return, [min(cloud[0,*]), max(cloud[0,*]), min(cloud[1,*]), max(cloud[1,*]), min(cloud[2,*]), max(cloud[2,*])]
end

function findboundsc, cloud
	a = findbounds(cloud)
	return, [a, a[0] + (a[1] - a[0]) / 2, a[2] + (a[3] - a[2]) / 2, a[4] + (a[5] - a[4]) / 2]
end

function findboundscabs, cloud
	a = findboundsc(cloud)
	return, [a, a[1] - a[0], a[3] - a[2], a[5] - a[4]]
end

function currentds, tid
	return, (( (_IDLitSys_GetSystem())->$
		getByIdentifier(tid))->$
			getByIdentifier("WINDOW/VIEW_1/VISUALIZATION LAYER"))->$
				GetCurrentDataSpace()
end

function drawcr, chm, seeds, objcrowns
	iimage, chm, identifier=t
	iplot, overplot=t, seeds, /scatter, sym_color=[255,0,0], use_default_colour=0
	(currentds(t))->add, objcrowns
	return, t
end

pro testtri, input
	bounds = findbounds(input)
;	plot, findgen(bounds[1] - bounds[0]), /nodata, xrange=[bounds[0], bounds[1]], yrange=[bounds[2], bounds[3]]
;	plots, input[0,*], input[1,*], psym=1

	triangulate, input[0,*], input[1,*], triangles, boundary

	s = size(triangles, /dimensions)
	num_triangles = s[1]
	Plot, findgen(bounds[1] - bounds[0]), /NoData, xrange=[bounds[0], bounds[1]], yrange=[bounds[2], bounds[3]]
	;PlotS, input[0,*], input[1,*], PSym=1
	FOR j=0,num_triangles-1 DO BEGIN
	   thisTriangle = [triangles[*,j], triangles[0,j]]
	   PlotS, input[0,thisTriangle], input[1,thisTriangle]
	ENDFOR

;	mesh_obj, 0, vert_list, poly_list, input
;
;	iplot, /scatter, input[0:1,*], identifier=tid
;	vis = obj_new('IDLitVisPolygon')
;	vis-> setproperty, data=input[0:1,*], polygons=poly_list
;	vis -> setproperty, linestyle=1, fill_background=0, noClose=0, shading=0
;
;	(currentds(tid))->add, vis
end

pro overplotAll, chm, objcrowns, tid
	iimage, chm, overplot=tid
	(currentds(tid))->add, objcrowns
end

;output is 2,seeds * numextents list of vertices, numextents vertices belong to each seed
function makeverts, seeds, extents, cycle
	numseeds = n_elements(seeds[0,*])
	numextents = n_elements(extents[0,*])
	verts = lonarr(numseeds * numextents,2)
	for i=0l, numseeds - 1,4 do begin
		verts[i:i+3,*] =  transpose(rebin(seeds[*,i], numextents, 2)) + (cycle[*,*] * transpose(rebin(extents[i,*], numextents, 2)))
	endfor
	return, verts
end

function buildpolyline, ncirc, nvertspc
	conn = ulonarr(ncirc * (nvertspc + 1))
	for i=0L, ncirc - 1 do begin
		conn[i * (nvertspc + 1):(i+1)*(nvertspc + 1)-1] = [nvertspc, indgen(nvertspc) + (i * nvertspc)]
	endfor
	return, conn
end

;nvertspe typically 3 or 4
;converts an z[n x m] array to a 1D array of form: [n, z[n,m],...]
function array_to_connectivity, input
	sz = size(input, /dimensions)
	conn = ulonarr((sz[0] + 1) * sz[1])
	for i=0l, sz[1] - 1 do $
		conn[i * (sz[0] + 1):(i+1)*(sz[0] + 1) - 1] = [sz[0], input[0,i], input[1,i], input[2,i]]
	return, conn

end


;ptr_list is an array of pointers to tree point arrays
function tree_alphas, points, clusters, num_clusters, valid_found=valid_found
	ptr_obj_poly_list = make_array(num_clusters, /ptr)
	ind = [0]
	validk = 0l
	for i=0l, num_clusters - 1 do begin
		ind = where(clusters eq i)
		if(n_elements(ind) ge 3) then ptr_obj_poly_list[validk++] = ptr_new(alphashape(points[*,ind],/auto_alpha))
	endfor
	if(arg_present(valid_found)) then valid_found=validk
	return, ptr_obj_poly_list
end

pro do_kmeans_plot, input, seeds, max_dist, max_iterations
	clusters = kmeans(input, seeds, max_dist, max_iterations, final_seeds=final_seeds)
	iplot, /scatter, identifier=tid, input
	iplot, /scatter, overplot=tid, seeds ,sym_color=[255,0,0]
	iplot, /scatter, overplot=tid, final_seeds, sym_color=[0,255,0]
	cds = currentds(tid)
	treeas = tree_alphas(input, clusters, n_elements(seeds[0,*]), valid_found=vf)
	for i=0L, vf -1 do cds->add, *(treeas[i])
end




function tracevalley, image, thresh
	out = image
	out[*] = 1

	minv = min(image, subs)
	xyind = array_indices(image, subs)

	out[where(image gt minv + thresh)] = 0

	return, out
end

function multiscale_maxima_old, image, wstart, maxwsize, increment=increment
	imout = image
	imout[*] = 0
	for i=wstart, maxwsize, 2 do begin
		imout += maxima(image, i)
	endfor
	return, imout
end



function is_maximum, vals, x
	return, ((x ge max(vals)) ? 1 : 0)
end

function readplots
	path = "D:\sjwelch\plotLoc2.txt"
	OPENR, lun, path, /GET_LUN

	lines = FILE_LINES(path)
	arrin = dblarr(2, lines)
	count = 0
	while(NOT(eof(lun))) do begin
		x = 0.0
		y = 0.0
		readf, lun, x, y
		arrin[0,count] = x
		arrin[1,count++] = y

	endwhile
	close, lun
	return, arrin

end

